Counting statistics of cotunneling electrons 
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We describe a method for calculating the counting statistics of electronic transport through 
nanoscale devices with both sequential and cotunneling contributions. The method is based upon a 
perturbative expansion of the von Neumann equation in Liouvillian space, with current cumulants 
calculated from the resulting nonMarkovian master equation without further approximation. As ap- 
plication, we consider transport through a single quantum dot and discuss the effects of cotunneling 
on noise and skewness, as well as the properties of various approximation schemes. 

PACS numbers: 73.23.Hk, 73.23.-b, 73.63.Kv, 42.50.Lc 



Cotunneling, the transfer of electrons via intermedi- 
ate "virtual" states, can be an important mechanism in 
the transport of electrons through quantum dots (QDs) 
P, In the Coulomb blockade regime, sequential tun- 
nelling processes are exponentially suppressed and, since 
it only suffers an algebraic suppression, cotunneling be- 
comes the dominant current-carrying mechanism. Exper- 
imental interest in cotunneling has remained high from 
the earliest experiments on metallic grains Q and large 
quantum dots through to more modern experiments 
on few-electron sing Ic- [I i, [3 and double- i] quan- 
tum dots. 

From a theoretical perspective, cotunneling refers to 
processes fourth-order in the coupling between the sys- 
tem and the leads. Such processes can be taken into 
account in a number of different ways, e.g. 0, [l^, [O, 
d, El, [ii, [H, [Hi. Most relevant here is the reabtime 
diagrammatic approach [H, [3, [H, Ell, in which higher- 
order tunnelling processes are incorporated into a master 
equation in a systematic fashion . This theory has been 
extensively developed and successfully applied to numer- 
ous transport problems: not just single QDs, but also 
double dots |l7|. quantum dot spin- valves carbon 
nanotubes [l9|, and QD-intcrferomcters [111 . Whilst such 
higher-order calculations have typically been restricted to 
the stationary current, and more recently, the shotnoise 
[H, [ill , much interest presently surrounds the full count- 
ing statistics (ECS) of the current, i.e. in current correla- 
tions beyond the second-order shotnoise [lil, [H, [2^ . The 
last few years has seen the advent of experiments capa- 
ble of detecting the pa ssag e of single electrons through 
QD systems [24. [25l [26l [27 . [28 1 and the experimental de- 
termination of ECS. Recently, measurements of the 15th 
cumulant were reported for a single quantum dot [2^. 

In this article we bring together several strands to in- 
vestigate the influence of cotunneling on ECS. We de- 
rive a fourth-order master equation for the reduced den- 
sity matrix of an arbitrary mesoscopic system using the 
Liouvillian-space perturbation theory of Ref. (29j and 
show how counting fields may be added in this formalism. 
Given that the resulting master equation is nouMarko- 
viaii, we employ the nonMarkovian formalism of Elindt et 
al. [301 to obtain expressions for the current cumulants. 

The calculation of the ECS for a nonMarkovian ME 



with cotunneling was described in Ref. [31|, but the ap- 
proach followed here is somewhat different. Braggio et al. 
[3l| calculate the cumulant generating function (CGE), 
and thus all the cumulants, rigorously to fourth order in 
the dot-lead coupling. Here, we calculate the Liouvillian 
rigorously to fourth-order and make no further approxi- 
mations in obtaining the cumulants. It is one of the aims 
of this paper to investigate the differences between the 
predictions of these two approaches. 

We use the transport through a single QD as our exam- 
ple system. We study first the single resonant level (SRL) 
model. Exact solutions exist for this model and this al- 
lows an evaluation of various approximation schemes. We 
also study the effects of interaction on transport with an 
Anderson model. Of the higher-order cumulants, we fo- 
cus here on the skewness as the first correlator beyond 
the shotnoise. We compare, both on a formal and a nu- 
merical level, with the work of Braggio et a l 13 il l, and 
with the shotnoise results of Thielmann et al [If 



I. TRANSPORT MODEL 

We begin by specifying the general transport set-up 
under consideration here. The total Hamiltonian H = 
Hj-cs -I- -ffs + ^ is composed of reservoir, system, and in- 
teraction parts. We write the system part in its diagonal 
basis Hs = -^al^}(^|j where \a) is a many-body sys- 
tem state of Na electrons. We consider a set of reservoirs, 
labelled with an index a that includes spin and any other 
relevant quantum numbers. We assume noninteracting 
reservoirs with Hamiltonian 



E 

k.a. 



(1) 



where Wka is the energy of the fcth mode in lead a, aua 
is a lead annihilation operator, and we have included the 
chemical potential of lead a, ^a, at this point for conve- 
nience. 

To ease book-keeping, we introduce a compact single 
index "1" to denote the triple of indices (^i, ki,ai). The 
first index = ± indicates whether a reservoir operator 
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is a creation or annihilation operator: 



Ci = 



flfclQl, Cl — - 

Leaving sums implicit, the reservoir Hamiltonian reads 

Hies = i^ka + l^a) 0.+kaO,-ka 



(2) 



{lui + ^i) aiaY(55i,+. 



(3) 



In equilibrium, the reservoir electrons are distributed ac- 
cording to the Fermi function 



1 



pu/kgT _|_ Y ' 



(4) 



which, since we include the chemical potential in Eq. ([T|) 
and assume a uniform temperature, is the same for all 
reservoirs. 

Single-electron tunnelling between system and reser- 
voirs is described by the Hamiltonian 



V 



tkamCll.ndm + ^kam^ln^kai 



/ V - <^^"'^-^k 
ka.m 



(5) 



where d„j is the annihilation operator for single-particle 
level m in the system, and tham is a tunnelling amplitude. 
We write this interaction as 



(6) 



with coefficients t+feam = ikam and t^k^m = *Lm' ^"^^ 
system operators in the many-body system basis 

= ^(aMm|a')(5 (Na - ^a' + 1) |a>(a'|, (7) 



and i-m = im - We have made explicit here the change in 
system charge induced by the operator. Although these 
operators only depend on ^i, we label them with the full 
"1" index for convenience: jxm = Jfim- 
At time t = we posit a separable total density matrix: 



pit = 0) = ps(0)p?e^„ 



(8) 



with the system in arbitrary state ps(io) and reservoirs 
in thermal equilibrium. 



II. LIOUVILLE-LAPLACE SPACE 

We now construct the elements required to perform 
our perturbation calculation in Liouville-Laplace space. 
In this section and the next, we follow Refs.[29, 32, 33|. 
The total density matrix evolves according to the von 
Neumann equation: 



p{t)^-i [H,p{t)]=Cp{t), 



(9) 



which defines the Liouvillian super-operator C = 
—i [H, • ]. This Liouvillian consists of three parts: 



— Cics + Cs + Cv 



(10) 



with £,.os = [-ffrcs,* ], = -i[Hs,»], and Cy 
—i [V, • ]. We write the interaction Liouvillian as 



p 



(11) 



where p = ± is a Keldysh index corresponding to the 
two parts of the commutator. Superoperators A and ,/ 
are defined through their actions on arbitrary operator 
O: For the reservoir, we have 



aiO, p = 



Oai, p = - ' 
and analogously for the system 

jlniO, p = + 



Ojlm, P 



(12) 



(13) 



By organising the elements of density matrices into 
vectors, superoperators, such as the Liouvillian, take the 
form of matrices. This is a particularly convenient repre- 
sentation for the system Liouvillian. We write a general 
system density matrix, ps = J2aia2 Pqiq2 |ai)(a2|, as the 
vector \ps)) = J2a Pa\4'a)), where the single index a cor- 
responds to the double (01,02) such that the "ket" \<f>a)) 
corresponds to |ai)(a2|. The action of the free system 
Liouvillian L5 on vector \4>a)) is 



Ls\<f>a)) = -i [Hs, |ai)(a2| 



(14) 



where = Ea^ — Ea^ defines the Bohr frequencies. The 
vectors | are therefore the right eigenvectors of Liou- 
villian Ls- The left eigenvectors, ((0a |, fulfil 



'{a\L, 



-iAa((0a| 



(15) 



and together with the right eigenvectors form a bi- 
orthonormal set: {{(j>a\4>a')) = Sa.a'- We have the com- 
pleteness relation in Liouvillian space 



(16) 



In general, it is important to make the distinction be- 
tween left and right eigenvectors because an arbitrary 
super-operator, in particular the effective system Liouvil- 
lian, will not be Hermitian, and the left and right eigen- 
vectors are therefore not adjoint. 



III. EFFECTIVE LIOUVILLIAN 



With the definition of the Laplace transform 



piz)= / dte-^'pit), 
Jo 



equation ([9]) yields the solution 

1 



-C 



piO). 



(17) 



(18) 



3 



Tracing Eq. (jlSp over reservoir degrees of freedom, results 
in an expression for the reduced density matrix of the 
system that wc write 



Ps{z) = 



1 



z-W{z) 



Ps{to) 



(19) 



where yy(z) is the nonMarkovian effective dot Liouvil- 
han. This we write as 



w(z)-/:s + s(z), 



(20) 



with Cs describing the free evohition of the system and 
^(z) the self-energy or "memory kernel" arising from cou- 
pling with the leads. 

In the pcrturbativc approach pursued here, the mem- 
ory kernel is calculated as the series S(z) = ^''"''(^) 
where n corresponds to the number of interaction Li- 
ouvillians Cv incorporated in that term. Tunnelling is 
governed by the rates 

r^^'^'iuj) = 27TY,tl^JlmA^~^kc.), (21) 
fei 

the diagonal elements of which are the familiar Fermi 
golden rule rates 



ki 



27ry^ |<l,„iP(5(tJ - LOka)- 



(22) 



In these terms, the expansion of ^(z) is seen as an expan- 
sion in small parameter T = T/ksT, such that S(")(z) 
is order (F)"/^. In the current work, we expand up to 
fourth order in the coupling Hamiltonian (second order 
in F), such that 



E(z) «E(2)(z) + E(4)(z). 



(23) 



The first term describes sequential tunnelling and the 
second cotunneling. 

Details of the calculation of the memory kernel terms 
are given in Appendix[Xl Assuming a constant tunnelling 
density of states F(cj) = F, the sequential term reads 

v(2)/ \ _ jP2 ZEl 7-Pi 

Z — lt,2\yi2 + M2j — J-S 

Xt2rn,tlmi72r\ (24) 

where 7^1^^ (^2^A^^)cq. is an equilibrium reservoir 
correlation function which evaluates as 



(25) 



Wc then obtain 



-piP2JZJMUa\JZ, 

xr^r^a'n^;6,Pi,m), (26) 

with the regularised integral 

/(2)(z = 0+-*e) = i/(pi(A,+eiM«, -e)) 



in terms of the function 

0(A) = Rc^^(l-+^, ^ 



2 2TTkBT 



In- 



D 



2TTkBT 



(27) 



with 5* the digamma function. 

The cotunneling term has two contributions: "direct" 
and "exchange", such that S(4)(z) = J:(^d\z)+J:<-'^^\z). 
The direct part is given by 

E(4^)(Z) ^ P4PlJLj0a))((0a|J|:,3l0a')) 

xUa'\JZJM)Ua"\JZ, 
X^ 1 ^2 

X-^fa'a" il,^2,Pl,P2, ^1,^2), (28) 

and the exchange 

E(4^)(Z) = -p4Pl44J</'a))((0a|4:,J0a')) 
Xi ^ ±2 

Xiaa'a" (^ ; , 6 , Pi , P2 , Ml , M2) ■ (29) 

These fourth-order integrals, and , are discussed 
in the Appendix. 

IV. FULL COUNTING STATISTICS 



The density matrix of Eq. (|T9)l is the Laplace- 
transform of the solution to the nonMarkovian master 
equation [sl] 



ps(t)= / dt'W{t-t')p{t'). 
Jo 



(30) 



In contrast to e.g. Ref. [29[, we make no further approx- 
imations at this point — in particular, we do not make 
Markov approximation — but rather seek to calculate 
the FCS of Eq. ([501) it stands. This requires that we 
introduce counting fields Xa in the appropriate places in 
the Liouvillian, a process which yields the "x-resolved" 
Liouvillian W{x,z) [H. 

In the current Liouville scheme, this can be achieved 
by replacing each contraction "f2i^^ in the memory kernel 
S(z) by the counting-ficld-dependent analogue 7fi^^(x)j 
which we define as 



7frnx) 



P2P^ 

721 cxp 



Pi -P2 



(31) 



with Sq = ±1 a factor given by the sign-convention for 
current flow in lead a. In the following we will only 
count in a single lead, for which we choose Sq = 1. To 
see that 721''^ (x) adds a counting field at the correct 
points, consider, for example, the case with —^2 = ^i = 
+. Then, for — p2 = Pi = +1 the contraction ^2?^ is 
proportional to trace of ci^^paka, a state with one more 
electron in lead a than p itself. On the other hand, for 
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P2 = Pi = +, we require the trace of a\^akaP, a state 
with the same number of electrons as p. The required 
counting- field factors are e**°-^° and e° = 1, respectively, 
as given by Eq. ([51]) . 

Once in possession of the "x-resolved" Liouvillian, the 
cumulant generating function (CGF) J-{x) ~ ^^t^oix) is 
obtained from the solution zq(x) of the equation 



where 



zo - Ao(x; ^o) = 0, 



(32) 



where Ao(x; zq) is the eigenvalue of W(x; z) that develops 
adiabatically from zero as x is increased from zero [35|. 
In the Markovian case, Aq is independent of z, and the 
CGF is simply T{x) = — iAo(x). The nonMarkovian case 
is less straightforward, however. We follow here the ap- 
proach of Ref . d^l , which uses Eq. ([5^ to derive expres- 
sions for the cumulants themselves, by-passing an explicit 
evaluation of the CGF itself. This approach can deliver 
the cumulants up to, in principle, arbitrary order (> 20 
in Ref. (36j ) and without any further approximation. 
We first define 

Jix, e) = W(x, z = 0+-ie)- W(x = 0, z = 0+), (33) 



along with the derivatives 



(34) 



and analogously for higher-orders. We define the left and 
right null-vectors of W(0, 0+) via 



W(0,0+)|Vo)) = ((^o|W(0,0+) =0, 



(35) 



which we assume to be unique. The vector iV'o)) corre- 
sponds to the stationary density matrix of the system, 
and multiplication with (Lijin \ , corresponds to taking the 
trace over system states [3a]- We define the stationary 
state "expectation value" ((•)) = ((^/lol • IV'o)), and the 
projectors V — |'!/'o))((V'o| and Q = 1 — V. Finally, we 
require the pseudo-inverse 



7^(e) = Q 



1 



(36) 



From Refs. [30|, l3j 
are 



ie + L{0,0+ -ie) 
, the first three current cumulants 



{{I')) 

an 



an 



(37) 
(38) 



3((/')) 



ii{{i^))rn{{j" - 2j'nj' - j"nj)) 
&i{{i^)),n{{J'n [njvj' + j'))) 

3i (((/'))™)' iJ' + ij'njnj)) 
{{{i'))mY {{J'nj + 2j'nj)). 



iiJ')) 

{{J" - 2j'nj')) 
iJ'" - ?.j"nj' 



(40) 
(41) 

- ij'nj")) 
j'n)j')), (42) 



(39) 



are the cumulants in the Markov approximation. In these 
expressions, it is understood that the pseudo-inverse is 
evaluated at e = 0. Although it is only practicable to 
explicitly write down the cumulants up to third order, the 
high-order cumulants can be obtained recursively [30|. 

Reference [3l| took a different approach to calculating 
the cumulants. There it was assumed that W is known 
to a given order in some small parameter, and the CGF 
is then calculated to the same order. For problems such 
as that considered here, this means that the CGF, and 
hence all the cumulants, are calculated rigorously up to 
order (F)^. This is to be contrasted with the above cumu- 
lants which, if expanded, have contributions at all orders 
in (F). Whilst the method of Ref. [3l[ may seem more 
consistent, there are two good reasons why the approach 
described here might be preferable. Firstly, from Ref. [Ill 
we know that in the infinite bias limit, the effective Liou- 
villian is given exactly by £0 plus the rate part of E^^-* . In 
order to recover the FCS correctly in this limit then, no 
further approximations should be made when calculating 
the cumulants [s^. Secondly, Ref. [2^ makes the point 
that, in certain circumstances, by treating incoming and 
outgoing processes unequally, a strict order-by-order ap- 
proach can lead to unphysical results for the current. In 
the next section we shall compare these two methods di- 
rectly. 

Before doing so, we note that comparing the forego- 
ing expressions for the current and the shotnoise with 
those of, e.g., Ref. [l^, we can identify the derivatives 
of J with respect to x and z with the various "current 
super-operator blocks" of the real-time diagrammatic ap- 
proach. For example, differentiating J once with respect 
to ix and setting x ~^ yields a super-operator like W, 
but with an additional forefactor ^i(pi — P2)l2. Under 
the summation, the p2 contribution cancels, leaving a 
forefactor ^ipi/2. This forefactor is the same as arises in 
replacing a single tunnel vertex with a current vertex in 
the self-energy, which is the recipe for obtaining the cur- 
rent super-operator block used to calculate the stationary 
current in the diagrammatic approach. Using a very dif- 
ferent approach, we thus reproduce the real-time current 
and shotnoise expressions. The advantage of the present 
method is that it is now easy to obtain to the higher cu- 
mulants, whereas with the diagrammatic approach, this 
requires some effort. 
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V. TRANSPORT THROUGH A SINGLE 
QUANTUM DOT 

We model the transport through a single Zeeman-split 
level with the Anderson Hamiltonian 14011 



(43) 



where is the energy of a spin-cr electron in the dot 
and U is the interaction energy. The reservoir and in- 
teraction Hamiltonians are as Eq. ^ and Eq. with 
index a including both lead L, R) and spin index. In 
the limit of large level-splitting we can address one and 
only one Zeeman level. We then recover the SRL model, 
the transport properties of which can be obtained exactly 
from scattering theory [4l|, [4^ . The SRL thus provides a 
useful benchmark against which to compare our approx- 
imate methods. 

We will discuss results obtained in several different ap- 
proximate schemes. We denote as "0(4)" the results ob- 
tained by calculating the cumulants as outlined in the 
previous section with the fourth-order effective Liouvil- 
lian containing both sequential and cotunneling terms. 
The second-order "0(2)" solution is obtained in the same 
way, but with sequential terms only. We also consider a 
scheme in which we expand the cumulants and truncate 
at fourth order. In this way we recover the ECS results 
of Braggio et al [3l| and the shotnoise results of Thiel- 
mann et al [l^. This approach we label as "0(4) trunc" . 
Finally, we compare with results in the Markovian ap- 
proximation, which we label with "Mark." . 

We calculated results with and without the level- 
renormalisation parts of the fourth-order self-energy (in- 
tegrals and /-^°, and /^^ and /^^ of the Appendix). 
Their contribution was found to be negligible in all cases 
studied here. In the results presented below, these parts 
of the self-energy have been neglected. This reduces 
considerably the computational effort involved since the 
double-principal-part integrals {I^'^ and I-^'^) must be 
evaluated numerically. 



A. Single resonant level 

The calculation of the first three cumulants in the scat- 
tering approach is discussed in Appendix [BJ In the infi- 
nite bias limit, we have 



{I) 

5(3) 



_tR. 



ill 



^ (Ft - 2FiFfl 



p2 ' 

eFlF^j - 2FiF|j 



or (/) = Fl/2, S = Tl/4: and S'-^^ = Fl/8 for symmetric 
couphng, Fl = F_R. 

Figure [1] shows the current and shotnoise as a func- 
tion of applied bias for different values of the coupling 
f = Fi/fcsT = Vn/kBT. A step occurs at eV « 2e 




20 40 60 

eV/kBT 



20 40 60 

eV/ksT 



FIG. 1: Stationary current I = {I)/kBT (left) and zero- 
frequency shotnoise S = S/ksT (right) as a function of ap- 
plied bias eV for the single resonant level model with level lo- 
cated at e = 20fc_BT, chemical potentials ^ll = — = e.V/2, 
and bandwidth Dj= 10^ ksT. Results are shown for three dif- 
ferent couplings: F = TL/kBT = Tn/kBT = ^, |, 1 and three 
calculational schemes: exact, full 4th-order (0(4)), and 2nd- 
order Markovian (0(2) Mark.). Whereas the 0(2) Marko- 
vian results show obvious deviations from the exact results 
for these couplings, the 0(4) solution gives good agreement 
except around the top of the shotnoise step for F = 1. 



(= AOksT here) as the level enters the transport window. 
For F < i agreement between the our 0(4) fourth-order 
calculation and the exact result is excellent across the 
whole bias range. For greater couplings, F > 1/2, de- 
viation from the exact solution is seen in the shotnoise 
around the top of the step which signals the start of the 
break down of our approach. The difference between the 
second-order Markovian and the exact solution is stark. 
In the CB regime {eV < 30 here) the sequential current 
is almost totally suppressed, but the cotunneling current 
is still considerable. The 0(2) Mark, solution also shows 
significant error in the high-bias {eV > 2e) regime, which 
arises largely from the Markovian approximation. 

Figure [2] plots the skewness which show a pronounced 
undulation at the onset. The 0(2) Markovian solution 
provides only the coarsest description of this behaviour, 
whereas it is reproduced by the 0(4) solution. For T < j, 
the quantitative agreement with the exact result is good. 
Nevertheless, for a given coupling, the error is larger for 
skewness than for the noise. This is not too surprising, 
since we expect higher-order correlators to be sensitive 
to higher-order processes, which are of course neglected 
here. The implication is that for a given coupling, only 
cumulants up to a certain order can reliably be calculated 
from any approximate effective Liouvillian [4^ . Figures 
[3] and |4] compare the 0(4) and 0(4) trunc results at 
F = 1/2. We choose this value to highlight the differences 
between the two solutions which, for smaller couplings, 
are negligible. Near the top of the step (Fig. [3^ and 
Fig. 2^), the 0(4) and 0(4) trunc solutions are notice- 
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eV/keT 
45 50 



40 60 

eV/ksT 



80 



FIG. 2: Zero-frequency skewness S^^^ = S'^' /keT of SRL 
with the same parameters as in Fig. [1] For a given couphng, 
agreement with the exact solution is worse than for shotnoise, 
but still good at low couplings (e.g. F = 1/4). 



ably different, with our 0(4) results significantly closer 
to the exact result. Deep in the CB regime (Fig. [3]d), 
the two approximate solutions differ once more, but this 
time the trunc solution isjiiore accurate. The difference 
is small (a difference in S* of 2 x 10~^ near zero bias); 
it plays, however, a disproportionate role in determining 
the Fano factors in the CB regime as Fig. [5J; and Fig. [3)d 
show. From the exact solution, we know that noise Fano 
factor diverges at low bias (fluctuation dissipation theo- 
rem), whereas the skewness Fano factor F^^^^ = S'^^^ / (I) 
tends to unity. In this regime, both Fano factors are re- 
produced better by trunc solution than by solution 0(4). 
It should be borne in mind that, in obtaining the Fano 
factors in this region, one is dividing one very small quan- 
tity by another and thus even small absolute errors can 
effect Fano factors quite dramatically. As a warning not 
to take the finer details of the Fano factor too seriously, 
we observe that the 0(2) Mark, solution in the CB regime 
reproduces the exact Fano factors better than any of the 
fourth-order results. This is deceptive since the individ- 
ual current and shotnoise obtained with this method arc 
vastly different from the exact results. 



From this comparison we obtain a degree of confidence 
in the cumulants of up to at least third order for F = 
1/4. Our 0(4) solution appears to perform better in the 
region where levels are crossing the chemical potentials, 
whereas the 0(4) trunc solution gives better results in 
the CB regime. In this regime, the transport properties 
are effectively Markovian, which can be demonstrated by 
comparing fourth-order solutions with and without the 
Markov approximation.. 
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eV/kBT 



40 60 

eV/kBT 



FIG. 3: Zero-frequency shotnoise (left panels) and Fano fac- 
tor F — S/ (/) (right panel) of SRL as a function of applied 
bias eV. The tunnel rate is fixed at F = 1/2 here, but other- 
wise the parameters are as Fig. [1] In addition to the approx- 
imation schemes discussed in Fig. [TJ results are also shown 
here from a rigorous expansion of the cumulants to 4th order 
(0(4) trunc). Around the current step (panel (a)), the full 
0(4) solution describes the behaviour better than 0(4) trunc. 
However, in the low bias. Coulomb blockade regime (panel 
(b)), it is the 0(4) trunc solution that matches the exact so- 
lution better. Panel (c) illustrates the dangers of considering 
the Fano factor alone: although 0(4) trunc reproduces the 
exact Fano factor extremely well in the Coulomb blockade 
regime, so does the 0(2) Markovian solution, which we know 
gives the current and shotnoise individually extremely poorly. 



B. Anderson Model 

The current, shotnoise and Fano factor of the Ander- 
son model with cotunneling were investigated in Ref. , 
and, as Fig. [5] shows, the present calculation broadly re- 
produces these results. The situation in which the lower 
dot level lies below the transport window is of particular 
interest. As observed in Ref. [ll], increasing the applied 
bias results in a large peak in the Fano factor around 
the point where the upper dot level enters the transport 
window. The peak exists in the sequential tunnel limit, 
but its height, width, and location are markedly altered 
by inelastic cotunneling processes. No peak occurs in the 
shotnoise itself; only in the Fano factor is this feature vis- 
ible. Figure [6] shows our results for the skewness in this 
situation. It is immediately clear that the skewness Fano 
factor also shows a peak, and that this is even more pro- 
nounced than that of the noise. Furthermore, the skew- 
ness itself exhibits a sharp peak, as inset Fig. [61d shows. 
As for the noise Fano factor, the presence of cotunneling 
significantly reduces the height and overall area of the 
peak. 

This superPoissonian behaviour indicates a signifi- 
cantly bunched electron flow, which can nicely be ex- 
plained with the dynamical channel blockade model of 
[45j . in which a single level (here, the lower) is but 
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FIG. 4: Zero-frequency skewness (left) and associated Fano 
factor f = S''3V(^) (right) of the SRL model as a function 
of applied bias eV. Same parameters and labels as Fig. [3] 
Once again, the 0(4) solution performs better in the onset 
region. 



weakly coupled to the collector. In the simple sequen- 
tial picture of Ref . [45| , the shotnoise and skewness Fano 
factors are predicted to be F2 = {1 + p)/{l — p) and 
^3 = (1 + Ap + p'^)/{l —p)"^, where 1/p is parameter cor- 
responding to the number of ways in which the dot can 
be filled. With p = 1/3 (corresponding to three ways of 
filling the dot: from the left and right into the lower level, 
and from the left only into the upper level), we obtain 
F2 = 2 and F3 — 11/2, which arc almost exactly the val- 
ues obtained by our sequential 0(2) results at the tops 
of the peaks. Cotunneling reduces the heights of peaks, 
and good agreement with the dynamical channel block- 
ade model can be obtained with the choice p = 0.272. 

Both noise and skewness figures show results for 0(4) 
and 0(4) trunc solutions. At high bias, these solutions 
agree closely with one another and both predict the same 
position and widths of the Fano factor peaks. However, 
the heights of the peaks are given differently in the two 
approaches, with the 0(4) trunc peaks being somewhat 
higher. Differences between the two solutions are most 
pronounced at low bias, when the Coulomb blockade is 
in effect. For example, the 0(4) solution predicts a small 
subPoissonan dip before the superPoissonian peak, which 
is absent in the 0(4) trunc results. From our studies of 
the SRL model, we expect the 0(4) trunc prediction to be 
more accurate in this regime, and this is borne out by the 
fact that the skewness Fano factor in the 0(4) calculation 
does not tend to unity with decreasing bias. Conversely, 
based on the SRL results, we expect the height of the 
Fano factor peaks to be better described by the 0(4) 
solution, i.e. the lower of the two values. 



VI. CONCLUSIONS 

We have described a method for calculating the count- 
ing statistics of an arbitrary mcsoscopic system that takes 




eV/ksT 

FIG. 5: Current (b), noise (c) and Fano factor (a) for the 
Anderson model as a function of applied bias eV. Parameters 
were chosen as in [l^: Fl = Fr = jksT, fiL = —f^R = ^^V, 
e| = -ISksT, ei = SfcsT, U = 40fcsT and bandwidth D = 
W'kBT. A prominent peak is observed in the Fano factor 
around a bias such that the transport window includes the 
upper level whilst the lower level is still included. Cotunneling 
reduces the size of the peak. 



into account both sequential tunnelling and cotunneling 
of electrons. This is achieved by performing a perturba- 
tive expansion of the von Neumann equation in Liouville- 
Laplace space and adding counting fields to the reservoir 
contractions. Current cumulants are then obtained with- 
out further approximation with the pseudo-inverse ap- 
proach. In principle, cumulants of arbitrary order can be 
calculated in this fashion. However, as we expect higher- 
order cumulants to be sensitive to higher-order tunnel 
processes, there will be an inevitable reduction in accu- 
racy as the order of the cumulant increases. Use of the 
pseudo-inverse means that this method is applicable to 
systems of large size, unlike methods which explicitly re- 
quire the eigenvalue Ao(Xj z) of the effective Liouvillian. 

We have studied here transport through a single quan- 
tum dot with both sequential and cotunneling contribu- 
tions. Of particular interest is the single-resonant level 
model as it provides an exact case against which we 
can compare. Good agreement with the exact results 
was found for both shotnoise and skewness up to ra- 
tio between coupling rate and reservoir temperature of 
r < 1/4. Moreover, we also compared our results with 
those obtained from a rigorous expansion of the cumu- 
lants up to fourth order in the tunnel coupling. From this 
comparison we conclude that for intermediate biases, and 
in particular when system levels lie near the chemical po- 
tentials of the leads, the approach described here gives 
better results. In the Coulomb blockade regime, however, 
the rigorous fourth-order approach is better, but here the 
dynamics are effectively Markovian. In general, the dif- 
ference between these two different fourth-order approx- 
imations is far less than that between fourth and sec- 
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FIG. 6: Skewness (b) and skewness Fano factor F^^'^ (a) for 
the Anderson model with parameters and labelling as Fig. (5] 
Not only the skewness Fano factor, but the skewness itself 
show a large peak as the top level enters transport window. 



ond order approximations, and decreases with decreasing 
systcm-rcservoir coupling. 

Future work includes the study of transport models 
with internal quantum degrees of freedom, such as the 
double quantum dot, to understand how the interplay of 
cotunneling and internal dynamics effects counting statis- 
tics. 



APPENDIX A: DERIVATION OF EFFECTIVE 
SYSTEM LIOUVILLIAN 

Following Refs. [1^, [H, [111, we start by defining the 
system operators 



gko 



(Al) 



such that the interaction Hamiltonian of Eq. ^ can be 
written V = Ciai.gi. Correspondingly, in Liouville space 
we have 



with A as before, and G defined via 



1 1 -0.gi 



9iO, p = 

V = 



(A2) 



(A3) 



The object is a dot-space superoperator with matrix 
elements [SJ 



5ss5s 



1, Ns- Ns' = even 
p, Ns - = odd 



(A4) 



where, N^ is the number of electrons in state s. Note 
that G{ = paHiraJl^. 



The reduced density matrix of the dot is given by trac- 
ing out the electron reservoirs 



ps{z) = TrR I p(z)| 



Trr 



1 



z-C 



p(0) . (A5) 



This we expand in powers of Cv to obtain 

ps(z) = TrR, {(17o(z) + ^lo{z)CvMz) + • • •) p(0)} (A6) 

with free propagator ^q{z) = [z — >Crcs ~ ^ ■ With 
substitution of Eq. (jA2[) . a typical term of the expansion 
Eq. (jXe)) reads 



\i=i 



(-0" [\{^iPi] TrR { M^)<tP-AI-GI- . . . 

TP^AP^GP,'noiz)p{to)Y (A7) 



. . cr'' 



Evaluating the action of the superoperators we obtain 
a factor and, as the G operators also contain cr, 

they evaluate at different positions in the chain as 



Our typical term then looks like 

/ all \ /odd \ 



(A8) 



xTrR { no{z)AP"GP" . . . f^o(^)p(to)} , (A9) 

and the next task is to separate dot and reservoir degrees 
of freedom. For this we can use the dot-reservoir su- 
peroperator commutation relation A^Gy = —pp'G^iAP. 
We will also need the following reasons: TrR^j-cs = , 
'CrcsPjSs = and A^Crcs = {Crcs - Xi) A^ with 



(AlO) 



The commutation of the dot-operators through the free 
propagators therefore changes the argument of the prop- 
agator: 



APQoiz) = noiz + xi)Al 



(All) 



Bringing all the A operators to the right of the G op- 
erators generates a factor which exactly cancels with the 
first product in Eq. (|A9|) . Our term becomes 

(odd \ 
Hp, j Ti-R { ns{z)GP-ns{z^.,)Gl-,' . . . 

...GP'nsiz,)G{'ns{z)APr ...A{'p{0)} (A12) 
with free dot propagator 

1 



^s{z) 



z - Cs 



(A13) 
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and Zm, = z + YJi^m+i I < m < n - 1. 

The the reservoir expectation values, 
Trrcs{ AP"...Af p^q} = (AP"...^f )eq, can be 
evaluated with the rules of Wick's theorem in Liouvillc 
space, which read [1^ 

• decompose (. . .)ros into pair-contractions 

• add minus sign for each interchange of A 

• omit factor (^11°''^'^ P'') arising from a super- 
operators 

• each pair contraction then contributes a factor 

(Af ) = = S^TPlfc.^ (-CiPi^i)- (A14) 

With these rules, our typical term becomes 
(-^)"^!s(^„)G'f,"^^s(^«-l)G^-^.. 

\ dccomps J 

where the last factor indicates a sum over all pair decom- 
positions with the relevant Wick sign (—1)^^. 

Comparison of this expression with Eq. (|19p and 
Eq. ([201) allows us to identify the self-energy as 



n Virrcd. / 

xGP"l]s(^«-i)G^"_Y ■ ■ • GP^'ns{zi)GP\ 
where the sum is over irreducible contractions only. 

1. Second order 

At second order, there is only one contraction, and we 
have 

E(2)(,) = cr ^ ^ . Gi^^ir 

z + X2 — Ls 
r,P2 ~-Pl/ai(-6Pl^l) ^pi 
Z + -I- //aj - L-s 



= -Gf\cj^,,)){{4>a\G\^ 



Pl/ai(-6Pl^l) 
Z + i^l{uJl + ^ia^) + i^a 



Re-expressing the G superoperators in terms of J super- 
operators and tunnel amplitudes, we have 



17712 ''1™ 



/ai(-6Pl'^l) 



' Z + iil{0Jl + Mai) + i^a ' 



With rates defined as in Eq. ([2T|) with the constant tun- 
nelling density of states approximation, Eq. (|A16p be- 
comes 

xrrJ:"4'H^;6,Pi,Mi), (A17) 



with integral 



/(2) 



27r 



Z + i^l{uJi + ^ai)+i^a' 



where the bandwidth D is assumed to be much larger 
than all other energy scales in the problem. We regularise 
the integrals by setting ^ — > 0"*" — ze with e wholly real. 
We arc therefore left to evaluate the integral 



1 



2tt 



faA~ilPl^l) 



D 



' iO+ -ie + i^i (wi + fiai ) + i^a 



Use of Dirac's identity, (iO"'" + x) ^ = P[l/x] — iTr6{x), 
allows us to write 

with (f> defined in Eq. ((27)) . correct to order 1/D. 



2. Fourth order 

At fourth order, there are two linked contractions, 
(41)(32) and (42)(31), we we label "D" for direct and 
"X" for exchange. The direct contribution reads 

xGrf}s(^i)Gf7ff^73l''^ (A18) 



This evaluates as Eq. ([28)) with the integral 

duji I doj' 



(2^)2 A3 - Ai 



/(Plt^l)/(P2W2) 



(A16) 



{iO+ + uji +LU2 - A2)(i0+ +UJI - A3) 
+ (Ai ^ A3) (A19) 

with Ai = ClMai + ^a" - e, A2 = 6Mai + 6Ma2 + A^' - £, 

A3 = CiMai + Aq — e. Note that here the (Ai ^ A3) 
symbol includes the pre-integral forefactor. 
The exchange term is 

= -Gfns{zs)Gfns{z2) 

xGl^ns{z,)G\^^iriir\ (A20) 
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which yields Eq. ([29)) with exchange integral 



1 



(27r)2 A2 - A3 - Ai 
1 



duJi I duj2f{piUJi)f{p2UJ2) 

1 



iO+ + — Ai iO+ + 0^2 — A3 
1 



zO+ + UJl + UJ2 — X2 



1 



iO+ + wi + W2 — — Aj 



(A21) 



with Ai ^ ClA^ai + Aa" - e, A2 = ClMai + 6Ma2 + Aq/ - e, 

and A3 = S,2^J^a2 + Aa - e. 

The analytic evaluation of these integrals is discussed 
in Ref. [2§|. Let us note here that Dirac's identity can be 
used to split these integrals up in three parts, e.g. = 
jDo _|_ jDi jD2 ^ where the number in the superscript 
corresponds to the number of delta functions appearing 
in the integrand. Integrals with a single delta function 
give the cotunncling rates and are the most important 
for the current scheme. These are evaluated as 



tDI 



P1P2 F(A2,A3) -F(A2,Ai) 



(2^)2 
, Pi 



As 



i^(A3 



-Ai 
F(Ai) 



{2ny 



A3 — Ai 



(A22) 



with 



F(A', A) = ( - b{X') [0(-A) - cf>{X' - A)] 



--0(A)+0(A'-A)/(A) 



(A23) 



F{X) = -f 0(A), and where b{x) 



.-IT _ l^ 



is the 

Bose-Einstein distribution. The single delta-function Ex- 
change integral gives 

jxi ^ PiP2 / f(A2,Ai)-f(Ai+A3,Ai) 
(27r)2 A2 - A3 - Ai 

F(A2,A3)-F(Ai+A3,A3; 



+ 



A2 A3 — Ai 



(A24) 



The remaining parts give rise to renormalisation of the 
system levels, and whereas I^^ and I^^ can easily be 
evaluated analytically, the double principal-part integrals 
must be performed numerically. In the models studies 
here, however, these fourth-order renormalisation parts 
were unimportant. 



bias are given by the well-known expressions 
^ ^ hj dET{E)[fL{E)-,fR{E)] 



S = 



^JdE{ T{E)[fL{l - h) + fail - Jr)] 

+T{E)il-T{E)){fL~fRf}, 



where T{E) is the transmission probability of the device 
and fx is the Fermi function of lead X. We have set here 
e = h = 1 and define the correlation functions in agree- 
ment with those of FCS. The corresponding expression 
for the skewness at finite temperature and bias is less well 
known. However, from the results for the symmetrised 
correlator of (46j (the appropriate quantity here), we have 

^^'^ = 4ym. = ^/ dE{3S,oo{E) - SoooiE)} (Bl) 
with 

S^oo = (l-T)2/i(l-/i)(l-2/i) 

+T(l-r)/i(l-/i)(l-2/^) 
Sooo = {l-TffL{l-fLKl-2fL)+T{l-TfaLR 

+T2(1 - T)aiiL + T^fRil - /i?)(l - 2fR), 



and 



axY = /x(l-/x)(l-2/y) + /x(l-/y)(l-2/x) 
+/y(l-/A-)(l-2/x). (B2) 

In the infinite bias limit, /l = 1 and f ji ~ 0, wc recover 



■'SYM. 



J dE{T{E){l-T{E)){l-2T{E))}^ 



which is the more familiar expression for skewness in the 
scattering approach. 

From Ref. (4^ we know that the transmission proba- 
bility of the single resonant level at energy e with partial 
widths Tl and Tr is 



T{E) 



Li R 



{E-eY + {TL+TRY/A- 



(B3) 
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